program tedistboundscicentered, rclass
	* version control: centered does the recentered bootstrap
	syntax varlist(min=4 numeric) [if] [in], ///
	[tevalues(numlist)] [y0values(numlist)] [strata(varlist)] [coverage(real .95)] [wd]
	tempname deltacdf cate estimates ses percentile cateci deltacdfci deltacdfwdci wdpercentile cdeltacdfci
	tempfile bsreps
	
	tokenize `varlist'
	local y "`1'"
	macro shift
	local d "`1'"
	macro shift
	local z "`1'"
	macro shift
	local s "`*'"
	macro shift
	
	* construct expression lists
	if "`tevalues'"=="" local tevalues 0
	local tind=0
	local dlist
	if "`wd'"!="" local wdlist
	foreach t in `tevalues' {
		local ++tind
		local dlist `dlist' dlb`tind'=r(dlb`tind') dub`tind'=r(dub`tind')
		if "`wd'"!="" local dlist `dlist' wdlb`tind'=r(wdlb`tind') wdub`tind'=r(wdub`tind')
	}
	local numts=`tind'
	
	mat `deltacdfci'=J(`numts',10,.)
	mat colnames `deltacdfci'= tevalue lb ub impercentile sigl sigu lbbsl lbbsu ubbsl ubbsu
	if "`wd'"!="" {
		mat `deltacdfwdci'=J(`numts',10,.)
		mat colnames `deltacdfwdci'= tevalue lb ub impercentile sigl sigu lbbsl lbbsu ubbsl ubbsu
	}
	if "`y0values'"!="" {
		local yind=0
		local ytind=0
		local clist
		local cdlist
		foreach yval in `y0values' {
			local ++yind
			local clist `clist' clb`yind'=r(clb`yind') cub`yind'=r(cub`yind')
			foreach t in `tevalues' {
				local ++ytind
				local cdlist `cdlist' cdlb`ytind'=r(cdlb`ytind') cdub`ytind'=r(cdub`ytind')			
			}
		}
		local numys=`yind'
		mat `cateci'=J(`numys',10,.)
		mat colnames `cateci'= y0value lb ub impercentile sigl sigu lbbsl lbbsu ubbsl ubbsu
		
		local numyts=`ytind'
		mat `cdeltacdfci'=J(`numyts',11,.)
		mat colnames `cdeltacdfci'= y0value tevalue lb ub impercentile sigl sigu lbbsl lbbsu ubbsl ubbsu
		
	}
	bootstrap `dlist' `clist' `cdlist',reps(500) saving("`bsreps'",replace) strata(`d' `z' `strata') : tedistboundsiter`wd' `varlist' `strata',tevalues(`tevalues') y0values(`y0values') 
	mat `estimates'=e(b)
	
	mat `ses'=e(se)
	*compute Imbens-Manski critical values
	local tind=0
	local col=0
	foreach t in `tevalues' {
		local ++tind
		local ++col
		local lb=`estimates'[1,`col'] 
		local lbse =   `ses'[1,`col']
		local ++col
		local ub=`estimates'[1,`col']
		local ubse =   `ses'[1,`col']
		if `lbse'==. local lbse=0
		if `ubse'==. local ubse=0
		
		if `lb'>`ub' {
			local ubnew = `lb'
			local lb=`ub'
			local ub=`ubnew'
			mat `estimates'[1,`col']=`ub'
			mat `estimates'[1,`col'-1]=`lb'
			mat `ses'[1,`col']=`lbse'
			mat `ses'[1,`col'-1]=`ubse'
		}
		if `lbse'>0 | `ubse'>0 {
			qui mata imbensmanski((`lb',`lbse',`ub',`ubse',`coverage'),"`percentile'")
			local upper=100*`percentile'
			local lower=100*(1-`percentile')
			preserve
			qui use "`bsreps'",clear
			local lbind=`col'-1
			local ubind=`col'
			qui centile dlb`tind',centile(`lower')
			local lbbsl=r(c_1)
			qui centile dlb`tind',centile(`upper')
			local lbbsu=r(c_1)
			qui centile dub`tind',centile(`lower')
			local ubbsl=r(c_1)
			qui centile dub`tind',centile(`upper')
			local ubbsu=r(c_1)
			restore
		}
		else {
			local lbbsl=`lb'
			local lbbsu=`lb'
			local ubbsl=`ub'
			local ubbsu=`ub'
			scalar `percentile'=`coverage'
		}
		mat `deltacdfci'[`tind',1]=(`t',`lb',`ub',`percentile',`lbse',`ubse',`lbbsl',`lbbsu',`ubbsl',`ubbsu')
		if "`wd'"!="" {
			local ++col
			local wdlb=`estimates'[1,`col'] 
			local wdlbse =   `ses'[1,`col']
			local ++col
			local wdub=`estimates'[1,`col']
			local wdubse =   `ses'[1,`col']
			if `wdlbse'==. local wdlbse=0
			if `wdubse'==. local wdubse=0
			if `lb'>`ub' {
				local wdubnew = `wdlb'
				local wdlb=`wdub'
				local wdub=`wdubnew'
				mat `estimates'[1,`col']=`wdub'
				mat `estimates'[1,`col'-1]=`wdlb'
				mat `ses'[1,`col']=`wdlbse'
				mat `ses'[1,`col'-1]=`wdubse'
			}
			if `wdlbse'>0 | `wdubse'>0 {
				qui mata imbensmanski((`wdlb',`wdlbse',`wdub',`wdubse',`coverage'),"`wdpercentile'")
				local wdupper=100*`wdpercentile'
				local wdlower=100*(1-`wdpercentile')
				preserve
				qui use "`bsreps'",clear
				local wdlbind=`col'-1
				local wdubind=`col'
				qui centile wdlb`tind',centile(`wdlower')
				local wdlbbsl=r(c_1)
				qui centile wdlb`tind',centile(`wdupper')
				local wdlbbsu=r(c_1)
				qui centile wdub`tind',centile(`wdlower')
				local wdubbsl=r(c_1)
				qui centile wdub`tind',centile(`wdupper')
				local wdubbsu=r(c_1)
				restore
			}
			else {
				local wdlbbsl=`wdlb'
				local wdlbbsu=`wdlb'
				local wdubbsl=`wdub'
				local wdubbsu=`wdub'
				scalar `wdpercentile'=`coverage'
			}
			mat `deltacdfwdci'[`tind',1]=(`t',`wdlb',`wdub',`wdpercentile',`wdlbse',`wdubse',`wdlbbsl',`wdlbbsu',`wdubbsl',`wdubbsu')
		}
	}

	if "`y0values'"!="" {
		local lastcdfcol=`col'
		local yind=0
		foreach yval in `y0values' {
			local ++yind
			local lb=`estimates'[1,2*`yind'-1+`lastcdfcol'] 
			local lbse =   `ses'[1,2*`yind'-1+`lastcdfcol']
			local ub=`estimates'[1,2*`yind'+`lastcdfcol']
			local ubse =   `ses'[1,2*`yind'+`lastcdfcol']
			if `lb'>`ub' {
				local ubnew = `lb'
				local lb=`ub'
				local ub=`ubnew'
				mat `estimates'[1,2*`yind'+`lastcdfcol']=`ub'
				mat `estimates'[1,2*`yind'+`lastcdfcol'-1]=`lb'
				mat `ses'[1,2*`yind'+`lastcdfcol']=`lbse'
				mat `ses'[1,2*`yind'+`lastcdfcol'-1]=`ubse'
			}
			qui mata imbensmanski((`lb',`lbse',`ub',`ubse',`coverage'),"`percentile'")
			local upper=100*`percentile'
			local lower=100*(1-`percentile')
			preserve
			qui use "`bsreps'",clear
			local lbind=2*`yind'-1+`lastcdfcol'
			local ubind=2*`yind'+`lastcdfcol'
			qui centile clb`yind' ,centile(`lower')
			local catelbbsl=r(c_1)
			qui centile clb`yind' ,centile(`upper')
			local catelbbsu=r(c_1)
			qui centile cub`yind',centile(`lower')
			local cateubbsl=r(c_1)
			qui centile cub`yind',centile(`upper')
			local cateubbsu=r(c_1)
			restore
			mat `cateci'[`yind',1]=(`yval',`lb',`ub',`percentile',`lbse',`ubse',`catelbbsl',`catelbbsu',`cateubbsl',`cateubbsu')
		}
		local lastcatecol=2*`yind'+`lastcdfcol'
		local ytind=0
		foreach yval in `y0values' {
			foreach t in `tevalues' {
				local ++ytind
				local lb=`estimates'[1,2*`ytind'-1+`lastcatecol'] 
				local lbse =   `ses'[1,2*`ytind'-1+`lastcatecol']
				local ub=`estimates'[1,2*`ytind'+`lastcatecol']
				local ubse =   `ses'[1,2*`ytind'+`lastcatecol']
				if `lbse'==. local lbse=0
				if `ubse'==. local ubse=0
				if `lb'>`ub' {
					local ubnew = `lb'
					local lb=`ub'
					local ub=`ubnew'
					mat `estimates'[1,2*`ytind'+`lastcatecol']=`ub'
					mat `estimates'[1,2*`ytind'+`lastcatecol'-1]=`lb'
					mat `ses'[1,2*`ytind'+`lastcatecol']=`lbse'
					mat `ses'[1,2*`ytind'+`lastcatecol'-1]=`ubse'
				}
				if `lbse'>0 | `ubse'>0 {
					qui mata imbensmanski((`lb',`lbse',`ub',`ubse',`coverage'),"`percentile'")
					local upper=100*`percentile'
					local lower=100*(1-`percentile')
					preserve
					qui use "`bsreps'",clear
					local lbind=2*`ytind'-1+`lastcatecol'
					local ubind=2*`ytind'+`lastcatecol'
					qui centile cdlb`ytind' ,centile(`lower')
					local cdeltacdflbbsl=r(c_1)
					qui centile cdlb`ytind' ,centile(`upper')
					local cdeltacdflbbsu=r(c_1)
					qui centile cdub`ytind',centile(`lower')
					local cdeltacdfubbsl=r(c_1)
					qui centile cdub`ytind',centile(`upper')
					local cdeltacdfubbsu=r(c_1)
					restore
				}
				else {
					local cdeltacdflbbsl=`lb'
					local cdeltacdflbbsu=`lb'
					local cdeltacdfubbsl=`ub'
					local cdeltacdfubbsu=`ub'
					scalar `percentile'=`coverage'
				}
				mat `cdeltacdfci'[`ytind',1]=(`yval',`t',`lb',`ub',`percentile',`lbse',`ubse',`cdeltacdflbbsl',`cdeltacdflbbsu',`cdeltacdfubbsl',`cdeltacdfubbsu')
			}
		}
	}
	
	
	display _newline "Bounds on Treatment Effect CDF"
	*display _newline as result "   TE value" _skip(5) "CI lb" _skip(5) "CDF lb" _skip(5) "CDF ub" _skip(5) "CI ub"
	*mat li `deltacdfci',noheader nonames
	mat li `deltacdfci',noheader
	return matrix deltacdfci = `deltacdfci'
	if "`wd'"!= "" {
		display _newline "Williamson-Downs Bounds on Treatment Effect CDF"
		*display _newline as result "   TE value" _skip(5) "CI lb" _skip(5) "CDF lb" _skip(5) "CDF ub" _skip(5) "CI ub"
		*mat li `deltacdfwdci',noheader nonames
		mat li `deltacdfwdci',noheader
		return matrix deltacdfwdci = `deltacdfwdci'
	}
	
	if "`y0values'"!="" {
		display _newline as result " Conditional ATE given `y'(0):"
		*display _newline as result "  Y(0) value" _skip(4) "CI lb" _skip(5) "CATE lb" _skip(5) "CATE ub" _skip(5) "CI ub"
		*mat li `cateci',noheader nonames
		mat li `cateci',noheader
		return matrix cateci=`cateci'
		
		display _newline as result " Conditional TE CDF given `y'(0):"
		*display _newline as result "  Y(0) value" _skip(4) "CI lb" _skip(5) "CATE lb" _skip(5) "CATE ub" _skip(5) "CI ub"
		*mat li `cateci',noheader nonames
		mat li `cdeltacdfci',noheader
		return matrix cdeltacdfci=`cdeltacdfci'
	}
end


* mata programs
mata:
	void imbensmanski(real vector params,string scalar impercentile)
	{ 
		lb=params[1]
		lbse=params[2]
		ub=params[3]
		ubse=params[4]
		alpha=params[5]
		
		delta=ub-lb
		maxse=max((lbse,ubse))
		
		S = solvenl_init()

        solvenl_init_evaluator(S, &imcondition())
		solvenl_init_argument(S, 1, delta)
		solvenl_init_argument(S, 2,alpha)
		solvenl_init_argument(S, 3,maxse)
		solvenl_init_type(S, "zero")
		solvenl_init_technique(S, "newton")
		solvenl_init_numeq(S, 1)
		solvenl_init_startingvals(S, 1.8)
		solvenl_init_iter_log(S, "on")
		x = solvenl_solve(S)
		st_numscalar(impercentile,min((max((normal(x),alpha)),(1+alpha)/2)))
	}
	
	void function imcondition(real scalar x, real scalar value, delta,alpha,maxse)
    {
        x
		value = normal(x+delta/maxse) - normal(-x)-alpha
    }
end


